#include<stdio.h>
#include<stdlib.h>
#include<math.h>
#include <omp.h>
#include <mpi.h>
#include <time.h>
#define NUM_THREADS 4
#define BILLION 1E9

double Max (double*, int);

/* matriz = > acesso pela regiao de memoria 
  ex: matriz[l][c] = l*ORDEM_MATRIX + c;
      matriz[1][1] = 1*3 + 1 = 4 <= v[4] (quinto elemento)
 */
void CriaIdentidade(double *MatrizAEstrela, double *MatrizA , double *VetorBEstrela, double *VetorB, double *VetorX_k , int J_ORDER)
{ 
    int i,j;
    #pragma omp parallel for num_threads(NUM_THREADS) 
    /* Cálculo da Matriz AEstrela */
    for(i=0; i<J_ORDER; i++)
    {
        #pragma omp parallel for num_threads(NUM_THREADS)
        for(j=0; j<J_ORDER; j++)
	{
            if(i != j)
                MatrizAEstrela[i*J_ORDER + j] = MatrizA[i*J_ORDER + j]/MatrizA[i*J_ORDER + i];
            else
                MatrizAEstrela[i*J_ORDER + j] = 0;
	    
	    /*Cálculo do Vetor BEstrela*/
	    VetorBEstrela[j] = VetorB[j]/MatrizA[j*J_ORDER + j];
	    
	    /*Cálculo do Vetor X*/
	    VetorX_k[j] = VetorBEstrela[j];
        }
    }
    return;
}

void proc_Calculate_x(double *a,double *b, double *x,double *x_k1,double *VetorDif, int J_ORDER,int block, int myrank)
{
	int i,j;
	double aux;
	aux = 0;
	
	/*Cálculo de X k+1*/
	#pragma omp parallel for num_threads(NUM_THREADS)
        for(i=0; i< block; i++)
	{
	    #pragma omp parallel for num_threads(NUM_THREADS)
	    /* Calcula a matriz a * vetor k */
            for(j=0; j< J_ORDER; j++)
	    {
                if( i + (myrank*block) != j ){
                    aux += a[i*J_ORDER + j]*x[j];
                }
            }
            /* Calcula K1 */
            x_k1[i] = b[i] - aux;

	    aux = 0;
        }
        
        for(i=0; i< block; i++)
	{
	    /*Cálculo da Diferenca entre os valores de k e k+1*/
	    VetorDif[i] = fabs( x_k1[i] - x[i + myrank*block]);
	}
        return;
}

void update_x( double *VetorX_k, double *VetorX_k1, int J_ORDER)
{
	int i;
	#pragma omp parallel for num_threads(NUM_THREADS)
	/*Atualizacao do Vetor de X k por X k+1*/
        for(i=0; i<J_ORDER; i++){
            VetorX_k[i] = VetorX_k1[i];
        } 
  return;
}

void copy(double *rec,double *send,int n)
{
  int i;
  for(i=0 ; i<n;i++)
      rec[i] = send[i];
}

int main (int argc, char *argv[])
{
    
    MPI_Init(&argc,&argv);

    int myrank;
    int proc;
    int block;
    int i,j, count= 1, J_ORDER, J_ROW_TEST, J_ITE_MAX;
    double Mr= 10, aux=0, J_ERROR, rowTest=0;
    double  *MatrizA,*MatrizAEstrela, *VetorB, *VetorBEstrela,*VetorX_k, *VetorX_k1, *VetorDif, *Buffer;
    double *MatrizA2;
    double *VetorB2;
    
    time_t startTime, endTime;
    
    MPI_Comm_size(MPI_COMM_WORLD, &proc);
    MPI_Comm_rank(MPI_COMM_WORLD, &myrank);    

    if(myrank==0)
    {    
	
	/*Leitura da ordem da matriz*/
	scanf("%d", &J_ORDER);
	if(J_ORDER % proc != 0)
	{
	  printf("Quantidade de processo invalido\n");
	  exit(-1);
	}
	/*Alocacao de espaco*/
	MatrizA = (double*) malloc(J_ORDER*J_ORDER*sizeof(double*));
	MatrizAEstrela = (double*) malloc(J_ORDER*J_ORDER*sizeof(double*));
	VetorB = (double*) malloc(J_ORDER*sizeof(double));
	VetorBEstrela = (double*) malloc(J_ORDER*sizeof(double));
	VetorX_k = (double*) malloc(J_ORDER*sizeof(double));
	VetorX_k1 = (double*) malloc(J_ORDER*sizeof(double));
	VetorDif = (double*) malloc(J_ORDER*sizeof(double));
	/* Buffer auxiliar para guardar ultimo vetor k */
	Buffer = (double*) malloc(J_ORDER*sizeof(double));
	/*Leitura da fila a ser avaliada*/
	scanf("%d", &J_ROW_TEST);

	/*Leitura do erro permitido*/
	scanf("%lf", &J_ERROR);

	/*Leitura do nro maximo de iteracoes permitido*/
	scanf("%d", &J_ITE_MAX);

	/*Leitura da matriz A*/
	for(i=0; i< J_ORDER; i++){
	    for(j=0; j<J_ORDER; j++)
		scanf("%lf", &MatrizA[i*J_ORDER + j]);
	}
	
	/*Leitura do vetor B*/
	for(i=0; i< J_ORDER; i++){
	    scanf("%lf", &VetorB[i]);
	}
	
	CriaIdentidade(MatrizAEstrela,MatrizA ,VetorBEstrela,VetorB,VetorX_k,J_ORDER);
    }     
    MPI_Bcast(&J_ORDER, 1, MPI_INT, 0, MPI_COMM_WORLD);
    MPI_Bcast(&J_ROW_TEST, 1, MPI_INT, 0, MPI_COMM_WORLD);
    MPI_Bcast(&J_ERROR, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
    MPI_Bcast(&J_ITE_MAX, 1, MPI_INT, 0, MPI_COMM_WORLD);
   
    block = J_ORDER / proc;
    
    MPI_Bcast(&block,1,MPI_INT,0,MPI_COMM_WORLD);
    
    MatrizA2 = (double*) malloc(block*J_ORDER*sizeof(double));
    VetorB2 = (double*) malloc(block*sizeof(double));
    
    if(myrank!=0)
    {
	VetorX_k = (double*) malloc(J_ORDER*sizeof(double));
	VetorX_k1 = (double*) malloc(block*sizeof(double));
	VetorDif = (double*) malloc(block*sizeof(double));
    }
    
    if(myrank == 0)
      MPI_Bcast(VetorX_k,J_ORDER,MPI_DOUBLE,0,MPI_COMM_WORLD);
    else
      MPI_Bcast(VetorX_k,J_ORDER,MPI_DOUBLE,0,MPI_COMM_WORLD); /* recebe vetor */

    MPI_Scatter(MatrizAEstrela, block*J_ORDER, MPI_DOUBLE, MatrizA2, block*J_ORDER, MPI_DOUBLE, 0,MPI_COMM_WORLD);
    MPI_Scatter(VetorBEstrela,block,MPI_DOUBLE,VetorB2,block,MPI_DOUBLE,0,MPI_COMM_WORLD);
    MPI_Barrier(MPI_COMM_WORLD);
    if(myrank == 0)
    {
      printf("Inicio tempo...\n\n");
      startTime = time(NULL);
    }
    do{
       if(myrank==0)
       {
	  copy(Buffer,VetorX_k,J_ORDER);
       }
       proc_Calculate_x(MatrizA2,VetorB2,VetorX_k,VetorX_k1,VetorDif,J_ORDER,block,myrank);
       MPI_Barrier(MPI_COMM_WORLD);
       MPI_Allgather(VetorX_k1,block,MPI_DOUBLE,VetorX_k,block,MPI_DOUBLE,MPI_COMM_WORLD);
    
       if(myrank == 0)
       {
	  /*Com todos dados juntos calcula Mr*/
	  for(i=0;i<J_ORDER;i++)
	  {
	     VetorDif[i] = fabs(VetorX_k[i]-Buffer[i]); /* Note que utiliza um Buffer antigo de VetorX_k*/
	  }
	  /* Note que VetorX_k ja esta atualizado portanto X_K = X_K+1*/
	  Mr = Max(VetorDif,J_ORDER)/Max(VetorX_k,J_ORDER);
       }
      /* Envia sinal para todos que terminou */ 
      MPI_Bcast(&Mr,1,MPI_DOUBLE,0,MPI_COMM_WORLD);
      count++;
      MPI_Barrier(MPI_COMM_WORLD);
    }while( count < J_ITE_MAX && Mr > J_ERROR);
    
    if(myrank == 0)
    {
      /*Row Test*/
      for(j=0; j<J_ORDER;j++)
      {
	rowTest += MatrizA[J_ROW_TEST*J_ORDER + j] * VetorX_k[j];
      }
      endTime = time(NULL);
      printf("Tempo de execução: %.1lf segundos\n", difftime (endTime,startTime));
      printf("Interations:%d\n",count);
      printf("RowTest: %d => [%lf] =? %lf\n\n\n",J_ROW_TEST,rowTest,VetorB[J_ROW_TEST]);
    }
    MPI_Finalize();
    return 0;
}

/*Funcao para calcular o valor máximo de um vetor*/
double Max (double Vetor[], int J_ORDER)
{
    int i;
    double max;

    max = fabs(Vetor[0]);

    for(i=1; i<J_ORDER; i++)
    {
        if(fabs(Vetor[i])>max)
            max = fabs(Vetor[i]);
    }
    return max;
}